library(enrichplot)
library(ggpubr)
library(ggrepel)
library(pheatmap)
library(msigdbr)
library(tidyverse)
library(clusterProfiler)
source('00_util_scripts/mod_bplot.R')

# wstlkomic computed logfc ------------
wst_slevhc <- read_csv('mission/FPP/xiangya_sle_ms/wstlkomic/SLE_NC_all_volcano.csv')

wst_sle_fc <- wst_slevhc |>
  select(c(1,`log2(foldchange)`,P_value_adjust)) |>
  dplyr::rename('avg_log2FC' = 'log2(foldchange)',
                'p_val_adj' = 'P_value_adjust') |>
  separate_wider_delim(...1, '_', names = c('UNIPROT','gene'), too_many = 'drop')

wst_sle_fc <- wst_sle_fc |>
  pull(UNIPROT) |>
  bitr(fromType = 'UNIPROT', toType = 'SYMBOL',
       OrgDb = 'org.Hs.eg.db') |>
  right_join(wst_sle_fc) |>
  as_tibble()

wst_sle_fc <- wst_sle_fc |>
  mutate(gene = case_when(is.na(SYMBOL) & gene == 'NA' ~ UNIPROT,
                          is.na(SYMBOL) ~ gene,
                          .default = SYMBOL))

wst_sle_fc |>
  filter(p_val_adj < .05)

wst_sle_fc |>
  plot_bill_volc(exp_group = 'SLE', fc_thres = log2(1.5))

wst_sle_fc |>
  filter(gene %in% kegg_mva)

# ora ----------
## Gene Ontology ----------
upora <- wst_sle_fc |>
  filter(avg_log2FC > 0 & p_val_adj < .05) |>
  pull(UNIPROT) |>
  enrichGO(OrgDb = 'org.Hs.eg.db',
           keyType = 'UNIPROT',
           ont = 'BP',
           readable = T)

upora@result |>
  head(10) |>
  plot_enrichment()

upora <- upora |>
  simplify()

## kegg ----------------
wstkegg <- wst_sle_fc$UNIPROT |>
  bitr_kegg(fromType = 'uniprot', toType = 'kegg', organism = 'hsa')

upora_kegg <- wst_sle_fc |>
  filter(avg_log2FC > 0 & p_val_adj < .05) |>
  pull(UNIPROT) |>
  enrichKEGG(organism = 'hsa',
           keyType = 'uniprot')

upora_kegg@result |>
  filter(category == 'Organismal Systems') |>
  filter(qvalue < .05) |>
  as_tibble() |>
  plot_enrichment() +
  labs(title = 'Upregulated KEGG pathways in SLE proteome vs NC',
       subtitle = 'Category: Organismal Systems')

upora_kegg@result |>
  filter(category == 'Cellular Processes') |>
  filter(qvalue < .05) |>
  as_tibble() |>
  plot_enrichment() +
  ggtitle('Upregulated KEGG pathways in SLE proteome')

leuko_path <- upora_kegg@result |>
  as_tibble() |>
  filter(str_detect(Description, 'Leukocyte')) |>
  pull(geneID) |>
  str_split_1('/') |>
  bitr(fromType = 'UNIPROT',
       toType = 'SYMBOL',
       OrgDb = 'org.Hs.eg.db')

## MKEGG ---------
## KEGG Module
upora_mkegg <- wst_sle_fc |>
  filter(avg_log2FC > 0 & p_val_adj < .05) |>
  pull(UNIPROT) |>
  enrichMKEGG(keyType = 'uniprot')

upora_mkegg@result |> slice_min(qvalue)

## DAVID -----------
## enrich tool developed by NIH, data collected from KEGG, GO, reactome...
## need RDAVIDWebService package installed (already removed from latest bioc release)
upora_david <- wst_sle_fc |>
  filter(avg_log2FC > 0 & p_val_adj < .05) |>
  pull(UNIPROT) |>
  enrichDAVID(idType = 'uniprot', species = 'human')

## Pathway Commons ---------
## include many sources: reactome(166), pid(36), pathbank(0), panther(3), netpath(11), kegg(0), inoh(1), humancyc(2), Detailed(279), All(279)
## only support human genes
upora_pathcom <- wst_sle_fc |>
  filter(avg_log2FC > 0 & p_val_adj < .05) |>
  pull(UNIPROT) |>
  enrichPC(keyType = 'uniprot', source = 'All')
166+36+3+11+1+2==219
upora_pathcom@result |> count(qvalue < .05)

## WikiPathways ------------
## support multi species & only entrez id
upora_wikip <- wst_sle_fc |>
  filter(avg_log2FC > 0 & p_val_adj < .05) |>
  pull(UNIPROT) |>
  bitr(fromType = 'UNIPROT', toType = 'ENTREZID',
       OrgDb = 'org.Hs.eg.db') |>
  pull(ENTREZID) |>
  enrichWP(organism = 'Homo sapiens')

upora_wikip@result |> as_tibble()

# fig: heatmap -------
sle_sample <- wst_slevhc |>
  colnames() |>
  str_subset('SLE|NC')

wst_anno <- sle_sample |>
  str_extract('SLE|NC') |>
  as_tibble()

wst_anno <- tibble(name = sle_sample,
                   group = str_extract(sle_sample, 'SLE|NC')) |>
  column_to_rownames('name')

wst_leuko <- wst_slevhc |>
  separate_wider_delim(...1, '_', names = c('UNIPROT','gene'), too_many = 'drop') |>
  filter(UNIPROT %in% leuko_path$UNIPROT)

wst_leuko |>
  select(2:16) |>
  column_to_rownames('gene') |>
  pheatmap(scale = 'row', main = 'Leukocyte transendothelial migration pathway',
           annotation_col = wst_anno, fontsize_row = 8)

## mva ----------
kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps') |>
  str_to_upper()

wst_slevhc |>
  separate_wider_delim(...1, '_', names = c('UNIPROT','gene'), too_many = 'drop') |>
  filter(gene %in% kegg_mva) |>
  select(2:16) |>
  column_to_rownames('gene') |>
  pheatmap(scale = 'row', main = 'MVA metabolism pathway',
           cluster_cols = F, annotation_col = wst_anno)

# GSEA ----------
wst_sig_fc <- wst_sle_fc |>
  filter(p_val_adj < .05)

wst_sig_entrez <- wst_sig_fc |>
  pull(UNIPROT) |>
  bitr_kegg(fromType = 'uniprot',
            toType = 'kegg',
            organism = 'hsa') |>
  left_join(wst_sig_fc, join_by(uniprot == UNIPROT)) |>
  as_tibble()

# some duplicate entrez id emerge from bitr
wst_sig_entrez <- wst_sig_entrez |>
  distinct(kegg, .keep_all = T)

wst_gse_list <- wst_sig_entrez$avg_log2FC |>
  set_names(wst_sig_entrez$kegg) |>
  sort(decreasing = T)

## GO ---------
sle_gse_go <- wst_gse_list |>
  gseGO(ont = 'BP', OrgDb = 'org.Hs.eg.db',pvalueCutoff = .1)

# 231, all significant
sle_gse_go@result |>
  as_tibble() |>
  filter(p.adjust < .05)

sle_gse_go |> dotplot()

sle_gse_go@result |>
  slice_max(NES, n = 10) |>
  ggplot(aes(Description, NES, fill = p.adjust)) +
  geom_col() +
  coord_flip()

sle_gse_gosim <- sle_gse_go |>
  clusterProfiler::simplify()

# 97 after simplify
sle_gse_gosim@result |>
  as_tibble() |> DT::datatable()

sle_gse_gosim@result |>
  as_tibble() |>
  arrange(desc(NES)) |>
  mutate(rank = seq_along(Description)) |>
  relocate(rank) |>
  filter(ID == 'GO:1901224')

### fig: gsea plot -----------
sle_gse_gosim |>
  gseaplot2('GO:1901224', color = 'red', base_size = 6,
            title = 'Positive regulation of non-\ncanonical NFkB signaling\nNES=2.01 FDR=0.018')

publish_pdf('gsea_Nfkb.pdf')

### fig: heatmap ------------
# extract go gene set from msigdb
nfkb_noncanon <- msigdbr(subcategory = 'GO:BP') |>
  filter(gs_exact_source == 'GO:1901224') |>
  pull(gene_symbol) |>
  unique()

nfkb_noncanon <- wst_sig_fc |>
  filter(SYMBOL %in% nfkb_noncanon & avg_log2FC > 0)

wst_nfkball <- wst_slevhc |>
  separate_wider_delim(...1, '_', names = c('UNIPROT','gene'), too_many = 'drop') |>
  filter(UNIPROT %in% nfkb_noncanon$UNIPROT)

wst_nfkball |>
  select(2:16) |>
  column_to_rownames('gene') |>
  pheatmap(scale = 'row', main = 'Positive regulation of non-\ncanonical NFkB signaling pathway',
           annotation_col = wst_anno, fontsize = 6, filename = 'nfkb_all_heatmap.pdf',
           height = 2.3, width = 2.8, treeheight_col = 10, treeheight_row = 10,
           show_colnames = F, cutree_cols = 2)

## interested signaling module ------ 
calmod_path <- c('CALM1','CALM2','PPP3CC','NFAT5','NFATC1','NFATC2')
pyk2_path <- c('PTK2B','KRAS','HRAS','NRAS','MAPK3','MAPK1')

sig_prot <- wst_sle_fc |>
  filter(SYMBOL %in% pyk2_path & avg_log2FC > 0) |>
  pull(UNIPROT)

sig_prot_ratio <- wst_slevhc |>
  separate_wider_delim(...1, '_', names = c('UNIPROT','gene'), too_many = 'drop') |>
  filter(UNIPROT %in% sig_prot) |>
  select(gene:NC23) |>
  pivot_longer(-1) |>
  mutate(group = str_extract(name, 'SLE|NC'))

### PYK2-RAS-ERK1/2 pathway boxplot --------
pyk2_pval <- wst_sle_fc |>
  filter(SYMBOL %in% pyk2_path & avg_log2FC > 0) |>
  mutate(y.position = 1.35, group1 = 'NC', group2 = 'SLE',
         p_val_adj = case_when(p_val_adj < .05 ~ signif(p_val_adj, 2) |> as.character(),
                               .default = 'NS'),
         gene = fct_relevel(gene, pyk2_path))

sig_prot_ratio |>
  mutate(gene = fct_relevel(gene, pyk2_path)) |>
  ggplot(aes(group, value, color = group)) +
  geom_boxplot() +
  facet_wrap(~gene, scales = 'free_y') +
  expand_limits(y = c(0,1.4)) +
  scale_color_manual(values = c('blue','red')) +
  stat_pvalue_manual(data = pyk2_pval, label = 'p_val_adj',size = 2) +
  labs(x = 'Group', y = 'Normalized abundance',
       title = 'Protein abundance of Pyk2-Ras-ERK1/2 axis modules') +
  theme_pubr() +
  theme_jpub

publish_pdf('pyk2_protein_boxplot.pdf',width = 60)

### PYK2 NRAS pathway boxplot --------
pyk2_pval2 <- pyk2_pval |>
  filter(p_val_adj != 'NS')

sig_prot_ratio |>
  mutate(gene = fct_relevel(gene, pyk2_path)) |>
  filter(gene %in% c('PTK2B','NRAS')) |>
  ggplot(aes(group, value, color = group)) +
  geom_boxplot() +
  facet_wrap(~gene, scales = 'free_y') +
  expand_limits(y = c(0,1.4)) +
  scale_color_manual(values = c('blue','red')) +
  stat_pvalue_manual(data = pyk2_pval2, label = 'p_val_adj',size = 2) +
  labs(x = 'Group', y = 'Normalized abundance',
       title = 'Protein abundance of Pyk2-Ras-ERK1/2 axis modules') +
  theme_pubr() +
  theme_jpub

publish_pdf('pyk2_nras_protein_boxplot.pdf',width = 60)

### jitter with mean line
sig_prot_ratio |>
  mutate(gene = fct_relevel(gene, pyk2_path)) |>
  ggplot(aes(group, value, color = group)) +
  geom_jitter(height = 0, width = .1) +
  stat_summary(geom = 'crossbar', fun = mean) +
  stat_summary(geom = 'errorbar', fun.data = mean_cl_boot, width = .3) +
  facet_wrap(~gene, scales = 'free_y') +
  theme_pubr() +
  scale_color_manual(values = c('blue','red'))

### MAPK cascade
sle_gse_go@result |> as_tibble() |> filter(ID == 'GO:0000165')

sle_gse_go |>
  gseaplot2('GO:0000165', color = 'red', base_size = 6,
            title = 'MAPK cascade\nFDR=0.74\nNES=0.936')

publish_pdf('gsea_mapk.pdf')

### ERK1 & ERK2
sle_gse_go@result |> as_tibble() |> filter(ID == 'GO:0070374')

sle_gse_go |>
  gseaplot2('GO:0070374', color = 'red', base_size = 6,
            title = 'Positive regulation of\nERK1 & ERK2 cascade\nFDR=0.29\nNES=1.39')

publish_pdf('gsea_erk12.pdf')

### calcineurin-NFAT signaling
sle_gse_gosim |>
  gseaplot2('GO:0033173', color = 'red', base_size = 6,
            title = 'calcineurin-NFAT signaling cascade')

publish_pdf('gsea_nfat.pdf')

sle_gse_gosim |>
  gseaplot2('GO:0070884', color = 'red', base_size = 6,
            title = 'Positive regulation of calcineurin-\nNFAT signaling cascade')

publish_pdf('gsea_nfatpos.pdf')

## KEGG --------
sle_gse_kegg <- wst_gse_list |>
  gseKEGG()

# 19, all significant
sle_gse_kegg@result |>
  as_tibble() |> DT::datatable()

## MKEGG --------
sle_gse_mkegg <- wst_gse_list |>
  gseMKEGG()
# no result

## Wikipath ------------
sle_gse_wp <- wst_gse_list |>
  gseWP(organism = 'Homo sapiens')

# 12
sle_gse_wp@result |> as_tibble() |> DT::datatable()

sle_gse_wp |>
  gseaplot2('WP3929', color = 'red', base_size = 6,
            title = 'Chemokine signaling pathway\nNES=1.91 FDR=0.022')

publish_pdf('gsea_chemokine.pdf')

### heatmap ----------
sle_chemo <- sle_gse_wp@result |>
  as_tibble() |>
  filter(ID == 'WP3929') |>
  select(core_enrichment) |>
  separate_longer_delim(1, '/') |>
  rename(kegg = core_enrichment) |>
  left_join(wst_sig_entrez)

wst_chemo <- wst_slevhc |>
  separate_wider_delim(...1, '_', names = c('UNIPROT','gene'), too_many = 'drop') |>
  filter(UNIPROT %in% sle_chemo$uniprot)

wst_chemo |>
  select(2:16) |>
  column_to_rownames('gene') |>
  pheatmap(scale = 'row', main = 'Chemokine signaling pathway',
           annotation_col = wst_anno, fontsize = 6, filename = 'chemokine_heatmap.pdf',
           height = 2.3, width = 2.8, treeheight_col = 10, treeheight_row = 10,
           show_colnames = F, cutree_cols = 2)

## PathCommons -----------
## only this func do not support entrez id
## seem to throw error

# modify function by remove error line
gsemy_pc <- function (geneList, source, keyType, ...) 
{
  keyType <- match.arg(keyType, c("hgnc", "uniprot"))
  source <- match.arg(source, clusterProfiler:::get_pc_source())
  pcdata <- clusterProfiler:::get_pc_data(source, keyType, output = "gson")
  res <- GSEA(geneList, gson = pcdata, ...)
  if (is.null(res)) 
    return(res)
  res@organism <- pcdata@species
  res@keytype <- keyType
  return(res)
}

sle_gse_pc <- wst_gse_list |>
  set_names(wst_sig_entrez$uniprot) |>
  gsemy_pc(source = 'All', keyType = 'uniprot')

# 71 for all, 63 for reactome
sle_gse_pc@result |> as_tibble() |> DT::datatable()

# zyx 1v1 data ==========
zyx <- read_csv('mission/zxy-prot-foldchange.csv', name_repair = 'universal')

zyx_fc2 <- zyx |>
  mutate(log2fc = log2(CKO.WT)) |>
  filter(abs(log2fc) > 1)

zyx_fc2_desc <- zyx_fc2 |>
  pull(log2fc) |>
  set_names(zyx_fc2$Accession) |>
  sort(decreasing = T)

## go all ------------
zyx_fc2_goall <- zyx_fc2_desc |>
  gseGO(ont = 'ALL',
        OrgDb = 'org.Mm.eg.db',
        keyType = 'UNIPROT')

zyx_fc2_goall@result |>
  as_tibble() |>
  write_csv('mission/zxy-ms-gsea-go.csv')

## KEGG ---------
zyx_fc2_entrez <- zyx_fc2 |>
  pull(Accession) |>
  bitr_kegg(fromType = 'uniprot',
            toType = 'kegg',
            organism = 'mmu') |>
  mutate(Accession = uniprot) |>
  left_join(zyx_fc2)

zyx_fc2_entrez_desc <- zyx_fc2_entrez |>
  pull(log2fc) |>
  set_names(zyx_fc2_entrez$kegg) |>
  sort(decreasing = T)

zyx_fc2_kegg <- zyx_fc2_entrez_desc |>
  gseKEGG(organism = 'mmu')

zyx_fc2_kegg@result |>
  as_tibble() |>
  write_csv('mission/zxy-ms-gsea-go.csv')

## path comm --------
zyx_fc2_reactome <- zyx_fc2_desc |> head()
  gsemy_pc(source = 'reactome', keyType = 'uniprot')
